﻿#include <memory>
#include "..\Storage\StorageAPI.h"
#include "..\Storage\Matrix.h"
namespace Preconditioning
{
        enum class PreconditionType
        {
				// LU - предобуславливание, U c единицами на диагонали. 
                LUfactor,
				// LDL - L с единицами на диагонали
				LDLfactor,
				// Диагональное предобуславливание
                DiagonalFactor,
				// QR - предобуславливание
                QRfactor,
				// LU(sq) - предобуславливание
                LUsqfactor,
				// LL - предобуславливание(неполное разложение Холесского)
				ILLfactor,
				// LU(*) - предобуславливание, у L на диагонали единицы
				LUwithunitL,
        };
       
		class Preconditioner : public Matrix
		{
		public:
			Preconditioner(){}
			virtual void MultPreconditionedDiagonalInverse(const Vector& to, Vector& res)const = 0;
		};
		template <typename BaseMatrix>
		class SimplePreconditioner :public Preconditioner
		{
		public:
			BaseMatrix base;
			SimplePreconditioner(int size = 0) :base(size){ flagL = false; flagU = false; }
			virtual void MultPreconditionedDiagonalInverse(const Vector& to, Vector& res)const override
			{
				res = to;
			}
			virtual void MultTransMatrixVector(const Vector& vec, Vector& res) const override
			{
				base.MultTransMatrixVector(vec, res);
			}
			virtual void multUpperTriagnle(bool withDiag, const Vector& vec, Vector& res) const override   //Aa?oiee
			{
				base.multUpperTriagnle(withDiag, vec, res);
			}
			virtual void multLowTriagnle(bool withDiag, const Vector& vec, Vector& res) const  override
			{
				base.multLowTriagnle(withDiag, vec, res);
			}
			virtual void multLInvToVector(bool withDiag, const Vector& in, Vector& out) const override     // (L [+D])^-1^ * vector
			{
				if (!flagL)
					withDiag = false;
				base.multLInvToVector(withDiag, in, out);
			}
			virtual void multUInvToVector(bool withDiag, const Vector& in, Vector& out) const override     // (U [+D])^-1^ * vector
			{
				if (!flagU)
					withDiag = false;
				base.multUInvToVector(withDiag, in, out);
			}
			virtual void getDiag(Vector& vec) const override //Iieo?eou aeaaiiaeu
			{
				throw new exception("bad function for preconditioner");
			}
			virtual int Size() const override
			{
				return base.Size();
			}
			virtual void Run(function<void(int, int, double)> fun) const
			{
				base.Run(fun);
			}
			virtual void multLTransInvToVector(bool withDiag, const Vector& in, Vector& out) const override// (L [+D])^T)^-1^ * vector, нужно для бисопряженных градиентов
			{
				base.multLTransInvToVector(withDiag, in, out);
			}

			virtual void multUTransInvToVector(bool withDiag, const Vector& in, Vector& out) const override// (L [+D])^T)^-1^ * vector, нужно для бисопряженных градиентов
			{
				base.multUTransInvToVector(withDiag, in, out);
			}

			friend static shared_ptr<Preconditioner> CreatePreconditioner(PreconditionType type, shared_ptr<Matrix> matr);
		private:
			// какую диагональ использовать
			bool flagL;
			bool flagU;
			// установить использование диагонали
			// (fL, fU - использовать диагональ нижнего и верхнего треугольника соответственно)
			virtual void setFlagDiag(bool fL, bool fU)
			{
				flagL = fL;
				flagU = fU;
			}
		};


		static shared_ptr<Preconditioner> CreatePreconditioner(PreconditionType type, shared_ptr<Matrix> matr)
		{
			typedef unsigned int uint;
			double zero = 1e-14;
			int Size = matr->Size();
			auto getElem([&](shared_ptr<Matrix> A, int i, int j)
			{
				double elem = 0;
				A->Run([&](int ii, int jj, double d)
				{
					if (i == ii && j == jj)
					{
						elem = d;
					}
				}
				);	
				return elem;
			});
			switch (type)
			{
//Автор - Жарков Константин
#pragma region LU - факторизация
			case Preconditioning::PreconditionType::LUfactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(Size));
				double sum = 0;
				double pr;
				double diag;
				Preconditioner->setFlagDiag(false, true);
				
				
				double elem;

				for (uint i = 0; i<Size;++i)
				for (uint j = 0; j<Size;++j)
				{
					elem = getElem(matr,i,j);
					if (fabs(elem) > zero)
					{
						if (i < j)
						{
							sum = 0;
							for (uint ik = 0; ik < j; ++ik)
								sum += getElem(Preconditioner, i, ik) * getElem(Preconditioner, ik, j);
							Preconditioner->base.SetElement(i, j, elem - sum);
						}
						if (i > j)
						{
							sum = 0;
							for (uint ik = 0; ik < j; ++ik)
								sum += getElem(Preconditioner, j, ik) * getElem(Preconditioner, ik, i);
							Preconditioner->base.SetElement(i, j, (elem - sum) / getElem(Preconditioner, j, j));
						}
						if (i == j && i != 0)
						{
							sum = 0;
							for (uint ik = 0; ik < i; ++ik)
								sum += getElem(Preconditioner, i, ik) * getElem(Preconditioner, ik, i);
							Preconditioner->base.SetElement(i, i, elem - sum);
						}
						if (i == 0 && j == 0)
							Preconditioner->base.SetElement(i, j, elem);
					}
				};
				return Preconditioner;
			}
#pragma endregion
#pragma region Диагональная факторизация
			case Preconditioning::PreconditionType::DiagonalFactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(Size));
				matr->Run([&](int i, int j, double d){if(i == j) Preconditioner->base.SetElement(i, i, sqrt(d)); });
				Preconditioner->setFlagDiag(true, true);
				return Preconditioner;
			}
#pragma endregion
#pragma region LDL факторизация
			case Preconditioning::PreconditionType::LDLfactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(matr->Size()));
				double sum = 0, diag, lower;
				Preconditioner->setFlagDiag(false, true);
				double el;
				//matr->Run([&](int i, int j, double el)
				for (uint i = 0; i < Size; ++i)
				{
					for (uint j = 0; j < Size; ++j)
					{
						el = getElem(matr, i, j);
						if (fabs(el) > zero)
						{
							if (j < i)
							{
								sum = 0;
								for (uint ik = 0; ik < j; ++ik)
									sum += getElem(Preconditioner, i, ik) * getElem(Preconditioner, ik, ik)  * getElem(Preconditioner, j, ik);// *getElem(Preconditioner, ik, ik);
								Preconditioner->base.SetElement(i, j, (el - sum) / (getElem(Preconditioner, j, j)));// *getElem(Preconditioner, j, j)));
							}
							if (i == j)
							{
								sum = 0;
								for (uint ik = 0; ik < i; ++ik)
								{
									lower = getElem(Preconditioner, i, ik);
									lower *= lower;
									diag = getElem(Preconditioner, ik, ik);
									//diag *= diag;
									sum += lower * diag;
								}
								Preconditioner->base.SetElement(i, i, (el - sum));
							}
						}
					}
				}
				return Preconditioner;
			}
#pragma endregion
			case Preconditioning::PreconditionType::QRfactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(matr->Size()));
				Preconditioner->setFlagDiag(false, true);
				Vector *uu = new Vector[Size];
				auto proj([&](Vector one, Vector two)
				{
					Vector vres;
					double res;
					res = one*two;
					res /= two*two;
					vres.resize(one.size());
					for (uint i = 0; i < one.size(); ++i)
						vres[i] *= res;
				});
				auto normalize([&](Vector u)
				{
					Vector res;
					double sum = 0;
					res.resize(u.size());
					for (uint i = 0; i < u.size(); ++i)
						sum += u[i];
					sum = sqrt(sum);
					for (uint i = 0; i < u.size(); ++i)
						res[i] = u[i] / sum;
					return res;
				});
				for (uint i = 0; i < Size; ++i)
				{

				}
				//return Preconditioner;
			}
				break;
// Автор - Морошкин Андрей			
#pragma region LU(sq) 
			case Preconditioning::PreconditionType::LUsqfactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(matr->Size()));
				Preconditioner->setFlagDiag(true, true);
				double sum;
				double el;
				//matr->Run([&](int i, int j, double el);
				for (uint i = 0; i < Size; ++i)
				{
					for (uint j = 0; j < Size; ++j)
					{
						el = getElem(matr, i, j);
						if (fabs(el) > zero)
						{
							if (i < j)
							{
								sum = 0;
								for (uint ik = 0; ik < j; ++ik)
									sum += getElem(Preconditioner, i, ik) * getElem(Preconditioner, ik, j);
								Preconditioner->base.SetElement(i, j, (el - sum) / getElem(Preconditioner, i, i));
							}
							if (i > j)
							{
								sum = 0;
								for (uint ik = 0; ik < j; ++ik)
									sum += getElem(Preconditioner, j, ik) * getElem(Preconditioner, ik, i);
								Preconditioner->base.SetElement(i, j, (el - sum) / getElem(Preconditioner, j, j));
							}
							if (i == j)
							{
								sum = 0;
								for (uint ik = 0; ik < i; ++ik)
									sum += getElem(Preconditioner, i, ik) * getElem(Preconditioner, ik, i);
								Preconditioner->base.SetElement(i, i, sqrt(el - sum));
							}
						}
					}
				}//);
				return Preconditioner;
			}
#pragma endregion

// Автор - Григоренко Семен
#pragma region LL - факторизация
			case Preconditioning::PreconditionType::ILLfactor:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(matr->Size()));
				int k;
				double sum, s, elem;
				Preconditioner->setFlagDiag(true, true);
				for (uint i = 0; i < Size; ++i)
				{
					for (uint j = 0; j < i+1; ++j)
					{
						elem = getElem(matr, i, j);
						if (fabs(elem) > zero)
						{
							if (i == j)
							{
								sum = 0;
								for (k = 0; k < i; k++)
								{
									s = getElem(Preconditioner, i, k);
									sum += s*s;
								}
								if ((elem - sum) >= 0)
									Preconditioner->base.SetElement(i, j, sqrt(elem - sum));
								else
									throw new exception("Sorry");
							}
							else
							{
								sum = 0;
								for (k = 0; k < j; k++)
								{
									s = getElem(Preconditioner, i, k);
									sum += s * getElem(Preconditioner, j, k);
								}
								s = getElem(Preconditioner, j, j);
								Preconditioner->base.SetElement(i, j, (elem - sum) / s);
							}
						}
					}
				}
				return Preconditioner;
			}
#pragma endregion

// Автор - Полищук Станислав
#pragma region LU(*) - факторизация, L - нижнетрегольная матрица с единицами на диагонали
			case Preconditioning::PreconditionType::LUwithunitL:
			{
				shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(Size));
				double sum = 0;
				double elem;
				Preconditioner->setFlagDiag(false, true);
				for (uint i = 0; i < Size; ++i)
				{
					for (uint j = 0; j < i; ++j)
					{
						elem = getElem(matr, i, j);
						if (fabs(elem) > zero)
						{
							sum = 0;
							for (uint k = 0; k < j; ++k)
							{
								sum += getElem(Preconditioner, k, j) * getElem(Preconditioner, i, k);
							}
							Preconditioner->base.SetElement(i, j, (elem - sum) / getElem(Preconditioner, j, j));
						}
					}
					for (uint j = i; j < Size; ++j)
					{
						elem = getElem(matr, i, j);
						if (fabs(elem) > zero)
						{
							sum = 0;
							for (uint k = 0; k < i; ++k)
								sum += getElem(Preconditioner, k, j) * getElem(Preconditioner, i, k);
							Preconditioner->base.SetElement(i, j, elem - sum);
						}
					}
				}
				return Preconditioner;
			}
#pragma endregion
			default:
				break;
			}
			shared_ptr<SimplePreconditioner<MappedRowSkyLineFormat>> Preconditioner(new SimplePreconditioner<MappedRowSkyLineFormat>(matr->Size()));
			Preconditioner->setFlagDiag(false, true);
			/*matr->Run([&](int i, int j, double elem)
			{
				if (i == j)
					Preconditioner->base.SetElement(i, j, 1);
			});*/
			for (uint i = 0; i < Size; ++i)
				Preconditioner->base.SetElement(i, i, (double)1);
			return Preconditioner;
			//return nullptr;
		}
}
